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Box  and  Muller  [1],  have  pointed  out  that  a  pair  of  independent 


standard  normal  random  variables  may  be  generated  in  the  form 
x  =  (-  2  In  u^)^cos  2tiu2,  y  =  (-  2  In  u^)^sin  2^^,  where  u^u^ 
are  independent  [0,1]  random  variables.  The  proof  readily  follows 
from  the  fact  that  (-  2  In  u^)^  is  distributed  as  p,  and  271U2  is 
distributed  as  6,  in  the  polar  representation  (p,9)  of  (x,y). 

The  method  is  attractive,  but  rather  slow,  requiring  sine,  cosine, 
logarithm,  and  square  root  subroutines.  This  note  shows  how  to  speed 
up  the  procedure  and  still  retain  its  principal  advantage  —  ease  of 
programming. 

One  line  of  improvement,  by  way  of  a  suggestion  of  von  Neumann  [2], 

2  2  -4  2  2  — i 

is  to  generate  sin  9  and  cos  9  in  the  form  u(u  +  v  )  ,v(u  +  v  )  , 

2  2 

where  (u,v)  are  uniformly  distributed  over  the  unit  circle:  u  +  v  <  1. 
Thi6  avoids  the  sine  and  cosine  subroutines.  The  pair  (u,v)  can  be 
produced,  with  efficiency  rc/4,  by  conditioning  a. pair  of  independent 
uniform  [-  1,1]  random  variables.  It  seems  then  that  we  may  avoid 
the  sine  and  cosine  subroutines  at  the  eepense  of  drawing  an  additional 
uniform  random  variable.  We  may,  however,  avoid  this  difficulty  by 
using  the  following  fact:  if]  u,v  are  uniformly  distributed  over  the 
unit  circle:  u  +  v  <  1,  then  u  +  v  is  uniform  [0,1]  and  is 
independent  of  u/v.  The  proof  is  elementary  and  omitted.  It 
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follows  that  u  +  v  is  independent  of  u(u  +  v  )  8  and  v(u  +  v  )  "f 


and  hence  that 


(1) 


x  =  u[-  2  ln(u2+v2)/(u2+v2)]^| 
y  =  v[-  2  ln(u2+v2)/(u2+v2)]^ 
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are  independent  standard  normal  random  variables,  since  [-  2  la(u  +  v  )]* 

2  2**i  2  2  —A- 

is  distributed  as  p,  and  u(u  +  v  )  ,v(u  +  v  )  8  are  distributed 

2  2 

as  sin  6,  cos  9,  and  are  independent  of  u  +  v  . 


We  then  suggest  this  procedure  for  generating  a  pair  x,y  of 
independent  standard  normal  random  variables:  generate  pairs  (u,v) 

of  independent  uniform  [-  1,1]  random  variables  until  one  satisfies 
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u  +  v  <  1 ,  then  form  x  and  y  according  to  relations  ( 1 ) .  This 
method  is  faster  than  the  direct  polar  coordinate  representation  and 
is  still  very  easy  to  program.  It  still  takes  about  5-6  times  as 
long  as  the  very  fastest  methods ,  but  may  well  serve  in  situations 
where  ease  of  programming  or  limited  storage  capacity  are  the  primary 
considerations. 


If  we  try.  to  generalize  the  above  methods  -  generate  a  point 
uniformly  over  the  surface  of  the  unit  n-sphere  and  project  the 
point  into  space  by  multiplying  by  a  chi-n  variate,  we  run  into  the 

r 

problem  of  producing  a  point  on  the  n-sphere.  The  obvious  method, 
generate  a  point  in  an  n-cube  and  reject  it  if  it  lies  outside  the 
inscribed  n-sphere,  has  efficiency  TTn^2/2nr(S^) ,  which  tends  rapidly 
to  zero.  If  we  do  not  specify  which  n-sphere  we  want,  however,  we  cah 
get  a  satisfactory  procedure  in  the  following  way:  generate  independent 
uniform  [-  1,1]  random  variables  vi>v2’  *  •  •  ,vn,vn+l  unt:1-1 

+”’+  vn  +  vn+r 


(2) 
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Then  the  point  ( v-j/S^ ,  , . . . ,  v ^/S  )  is  uniformly  distributed  over 

the  surface  of  the  unit  n-sphere,  since  the  density  function  of 

2  2 

(v^,...,vn),  given  condition  (2),  is  a  multiple  of  1  -  fl  -  (z^  +  •••  +  zr)  ] 
and  is  radially  symmetric . 


We  then  offer  this  method  for  generating  a  random  number  of 

independent,  standard  normal  random  variables  x^ ,x2» . . . ,xnt  Generate 

Independent  uniform  [-1,1]  random  variables  v^, . . . , vn>vn+1  until 

=  v2  +  •••  +  v2  <  1  <  v2  +•••+  v2  +  v2  ,  then  put  x.  =  v.  [r/s2]^, 
n  1  n  —  1  n  n+1  -  —  i  iL  v  nJ 

where  R  has  the  chl-square-n  distribution  and  is  independent 

of  the  v's.  If  n  =  2m,  we  may  put  R  =  -  2  ln(u^u2*  •  »um) 

2 

and  if  n  =  2m  +  1,  R  =  -  2  ln(u^u2*  •  ‘u^)  +  y  ,  where  y  is  normal. 

2 

We  may  produce  y  either  as  the  square  of  one  of  the  normal 
variates  previously  generated,  or  else  in  pairs  in  the  form 
w2  ln(w~  +  v2)/(w2  +  v2)  ,  w2  ln(w2  +  w2)/(v2  +  w2^  ’  where 
,w2  are  uniform  [0,1],  conditioned  by  w^  +  w,,  <  1. 


The  expected  number  of  x's  produced  by  this  method  is 
about  3*5. 
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